function [rmin,rmax] = approximateBounds(restr,irs,phi,opt)
% Approximate the bounds of the identified set using Algorithm 2.
% Inputs:
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

vma = irs.vmaGK;
H = opt.H;
ivar = opt.ivar;
jshock = opt.jshock;
KTilde = opt.KTilde;

% Draw from space of orthonormal matrices satisfying the equality and sign
% restrictions KTilde times.
opt.Qdraws = KTilde;
Q0 = drawQ0(restr,phi,opt);

etaDraw = zeros(H+1,length(ivar),KTilde); % Storage

for hh = 1:H+1 % For each horizon

    for ii = 1:length(ivar) % For each variable of interest

        cc = vma(ivar(ii),:,hh);
        
        for kk = 1:KTilde % For each draw of Q
            
            etaDraw(hh,ii,kk) = cc*Q0(:,jshock,kk); % Compute IRF
            
        end

    end

end

% For each horizon and variable, compute minimum and maximum IRF over draws
% of Q.
rmin = min(etaDraw,[],3);
rmax = max(etaDraw,[],3);

end

